{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# LU with Partial Pivoting"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 131,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import numpy.linalg as la\n",
        "\n",
        "np.set_printoptions(precision=3, suppress=True)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Set-up"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's grab a (admittedly well-chosen) sample matrix `A`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 132,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[  0.,   4.,  19.,  -7.],\n",
              "       [ -1.,  -2., -10.,  -0.],\n",
              "       [  1.,  17.,   1.,  -4.],\n",
              "       [ -5.,  -8.,  -6.,  -2.]])"
            ]
          },
          "execution_count": 132,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "n = 4\n",
        "\n",
        "np.random.seed(235)\n",
        "A = np.round(5*np.random.randn(n, n))\n",
        "A[0,0] = 0\n",
        "A[2,1] = 17\n",
        "A[0,2] = 19\n",
        "A"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "----------------\n",
        "## Permutation matrices\n",
        "\n",
        "Now define a function `row_swap_mat(i, j)` that returns a permutation matrix that swaps row i and j:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 133,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def row_swap_mat(i, j):\n",
        "    P = np.eye(n)\n",
        "    P[i] = 0\n",
        "    P[j] = 0\n",
        "    P[i, j] = 1\n",
        "    P[j, i] = 1\n",
        "    return P"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "What do these matrices look like?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 134,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.,  1.,  0.,  0.],\n",
              "       [ 1.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  1.,  0.],\n",
              "       [ 0.,  0.,  0.,  1.]])"
            ]
          },
          "execution_count": 134,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "row_swap_mat(0,1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Do they work?"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 135,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -1.,  -2., -10.,   0.],\n",
              "       [  0.,   4.,  19.,  -7.],\n",
              "       [  1.,  17.,   1.,  -4.],\n",
              "       [ -5.,  -8.,  -6.,  -2.]])"
            ]
          },
          "execution_count": 135,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "row_swap_mat(0,1).dot(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "--------------"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Part I"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "`U` is the copy of `A` that we'll modify:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 136,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "U = A.copy()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## First column"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Create P1 to swap up the right row: "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 137,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -5.,  -8.,  -6.,  -2.],\n",
              "       [ -1.,  -2., -10.,   0.],\n",
              "       [  1.,  17.,   1.,  -4.],\n",
              "       [  0.,   4.,  19.,  -7.]])"
            ]
          },
          "execution_count": 137,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "P1 = row_swap_mat(0, 3)\n",
        "U = P1.dot(U)\n",
        "U"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 138,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1. ,  0. ,  0. ,  0. ],\n",
              "       [-0.2,  1. ,  0. ,  0. ],\n",
              "       [ 0.2,  0. ,  1. ,  0. ],\n",
              "       [ 0. ,  0. ,  0. ,  1. ]])"
            ]
          },
          "execution_count": 138,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "M1 = np.eye(n)\n",
        "M1[1,0] = -U[1,0]/U[0,0]\n",
        "M1[2,0] = -U[2,0]/U[0,0]\n",
        "M1"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 139,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -5. ,  -8. ,  -6. ,  -2. ],\n",
              "       [  0. ,  -0.4,  -8.8,   0.4],\n",
              "       [  0. ,  15.4,  -0.2,  -4.4],\n",
              "       [  0. ,   4. ,  19. ,  -7. ]])"
            ]
          },
          "execution_count": 139,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U = M1.dot(U)\n",
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Second column\n",
        "\n",
        "Create `P2` to swap up the right row:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 140,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -5. ,  -8. ,  -6. ,  -2. ],\n",
              "       [  0. ,  15.4,  -0.2,  -4.4],\n",
              "       [  0. ,  -0.4,  -8.8,   0.4],\n",
              "       [  0. ,   4. ,  19. ,  -7. ]])"
            ]
          },
          "execution_count": 140,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "P2 = row_swap_mat(2,1)\n",
        "U = P2.dot(U)\n",
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make the second-column elimination matrix `M2`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 141,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.   ,  0.   ,  0.   ,  0.   ],\n",
              "       [ 0.   ,  1.   ,  0.   ,  0.   ],\n",
              "       [ 0.   ,  0.026,  1.   ,  0.   ],\n",
              "       [ 0.   , -0.26 ,  0.   ,  1.   ]])"
            ]
          },
          "execution_count": 141,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "M2 = np.eye(n)\n",
        "M2[2,1] = -U[2,1]/U[1,1]\n",
        "M2[3,1] = -U[3,1]/U[1,1]\n",
        "M2"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 142,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -5.   ,  -8.   ,  -6.   ,  -2.   ],\n",
              "       [  0.   ,  15.4  ,  -0.2  ,  -4.4  ],\n",
              "       [  0.   ,   0.   ,  -8.805,   0.286],\n",
              "       [  0.   ,   0.   ,  19.052,  -5.857]])"
            ]
          },
          "execution_count": 142,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U = M2.dot(U)\n",
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Third column\n",
        "\n",
        "Create `P3` to swap up the right entry:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 143,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -5.   ,  -8.   ,  -6.   ,  -2.   ],\n",
              "       [  0.   ,  15.4  ,  -0.2  ,  -4.4  ],\n",
              "       [  0.   ,   0.   ,  19.052,  -5.857],\n",
              "       [  0.   ,   0.   ,  -8.805,   0.286]])"
            ]
          },
          "execution_count": 143,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "P3 = row_swap_mat(3, 2)\n",
        "U = P3.dot(U)\n",
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Make the third-column elimination matrix `M3`:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 144,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.   ,  0.   ,  0.   ,  0.   ],\n",
              "       [ 0.   ,  1.   ,  0.   ,  0.   ],\n",
              "       [ 0.   ,  0.   ,  1.   ,  0.   ],\n",
              "       [ 0.   ,  0.   ,  0.462,  1.   ]])"
            ]
          },
          "execution_count": 144,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "M3 = np.eye(n)\n",
        "M3[3,2] = -U[3,2]/U[2,2]\n",
        "M3"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 145,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -5.   ,  -8.   ,  -6.   ,  -2.   ],\n",
              "       [  0.   ,  15.4  ,  -0.2  ,  -4.4  ],\n",
              "       [  0.   ,   0.   ,  19.052,  -5.857],\n",
              "       [  0.   ,   0.   ,   0.   ,  -2.421]])"
            ]
          },
          "execution_count": 145,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "U = M3.dot(U)\n",
        "U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Wrap-up"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "So we've built $M3P_3M_2P_2M_1P_1A=U$."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 150,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ -5.   ,  -8.   ,  -6.   ,  -2.   ],\n",
              "       [  0.   ,  15.4  ,  -0.2  ,  -4.4  ],\n",
              "       [  0.   ,   0.   ,  19.052,  -5.857],\n",
              "       [  0.   ,   0.   ,   0.   ,  -2.421]])"
            ]
          },
          "execution_count": 150,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "M3.dot(P3).dot(M2).dot(P2).dot(M1).dot(P1).dot(A)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "---------------------\n",
        "That left factor is anything but lower triangular:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 151,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.   ,  0.   ,  0.   ,  1.   ],\n",
              "       [ 0.   ,  0.   ,  1.   ,  0.2  ],\n",
              "       [ 1.   ,  0.   , -0.26 , -0.052],\n",
              "       [ 0.462,  1.   , -0.094, -0.219]])"
            ]
          },
          "execution_count": 151,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "M3.dot(P3).dot(M2).dot(P2).dot(M1).dot(P1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "# Part II"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Now try the reordering trick:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 160,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "L3 = M3\n",
        "L2 = P3.dot(M2).dot(la.inv(P3))\n",
        "L1 = P3.dot(P2).dot(M1).dot(la.inv(P2)).dot(la.inv(P3))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 155,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.   ,  0.   ,  0.   ,  1.   ],\n",
              "       [ 0.   ,  0.026,  1.   ,  0.   ],\n",
              "       [ 1.   , -0.266, -0.26 ,  0.   ],\n",
              "       [ 0.462,  0.878, -0.094,  0.   ]])"
            ]
          },
          "execution_count": 155,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "L3.dot(L2).dot(L1).dot(P3).dot(P2).dot(P1)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "--------------\n",
        "We were promised that all of the `L`*n* are still lower-triangular:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 168,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "[[ 1.   0.   0.   0. ]\n",
            " [ 0.2  1.   0.   0. ]\n",
            " [ 0.   0.   1.   0. ]\n",
            " [-0.2  0.   0.   1. ]]\n",
            "[[ 1.     0.     0.     0.   ]\n",
            " [ 0.     1.     0.     0.   ]\n",
            " [ 0.    -0.26   1.     0.   ]\n",
            " [ 0.     0.026  0.     1.   ]]\n",
            "[[ 1.     0.     0.     0.   ]\n",
            " [ 0.     1.     0.     0.   ]\n",
            " [ 0.     0.     1.     0.   ]\n",
            " [ 0.     0.     0.462  1.   ]]\n"
          ]
        }
      ],
      "source": [
        "print(L1)\n",
        "print(L2)\n",
        "print(L3)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "So their product is, too:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 172,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.   ,  0.   ,  0.   ,  0.   ],\n",
              "       [ 0.2  ,  1.   ,  0.   ,  0.   ],\n",
              "       [-0.052, -0.26 ,  1.   ,  0.   ],\n",
              "       [-0.219, -0.094,  0.462,  1.   ]])"
            ]
          },
          "execution_count": 172,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Ltemp = L3.dot(L2).dot(L1)\n",
        "Ltemp"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "----\n",
        "`P` is still a permutation matrix (but a more complicated one):"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 174,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.,  0.,  0.,  1.],\n",
              "       [ 0.,  0.,  1.,  0.],\n",
              "       [ 1.,  0.,  0.,  0.],\n",
              "       [ 0.,  1.,  0.,  0.]])"
            ]
          },
          "execution_count": 174,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "P = P3.dot(P2).dot(P1)\n",
        "P"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "-----------------\n",
        "So to sum up, we've made:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 175,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.]])"
            ]
          },
          "execution_count": 175,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "Ltemp.dot(P).dot(A) - U"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "--------------\n",
        "Multiply from the left by `Ltemp`${}^{-1}$, which is *also* lower triangular:"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 179,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 1.   ,  0.   ,  0.   ,  0.   ],\n",
              "       [-0.2  ,  1.   ,  0.   ,  0.   ],\n",
              "       [ 0.   ,  0.26 ,  1.   ,  0.   ],\n",
              "       [ 0.2  , -0.026, -0.462,  1.   ]])"
            ]
          },
          "execution_count": 179,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "L = la.inv(Ltemp)\n",
        "L"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 180,
      "metadata": {
        "collapsed": false
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([[ 0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.],\n",
              "       [ 0.,  0.,  0.,  0.],\n",
              "       [ 0., -0.,  0.,  0.]])"
            ]
          },
          "execution_count": 180,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "P.dot(A) - L.dot(U)"
      ]
    }
  ],
  "metadata": {},
  "nbformat": 4,
  "nbformat_minor": 0
}